
rm(list = ls())
## with mxmaps package https://www.diegovalle.net/mxmaps/

# if (!require("devtools")) {
#   install.packages("devtools")
# }
# devtools::install_github("diegovalle/mxmaps")

library(tidyverse)
library(lubridate)
library(forcats)
library(tidyr)
library("mxmaps")
library("ggrepel")
data("df_mxmunicipio")

here::here()

#### Prepare data
### With data on journalists murdered
journ = readxl::read_xlsx('Master Dataset_Journalists Killed.xlsx')
journ$year = year(journ$`Fecha (month-day-year)`)

murders_by_mun <- journ %>%
  # include leading zeros in variable munic_id
  dplyr::mutate(munic_id = str_pad(mun_id, 3, pad = "0"),
                st_id = str_pad(state_id, 2, pad = "0")) %>%
  # concatenate state & munic IDs in a in a single variable
  unite(CVEGEO, c("st_id", "munic_id"), sep = "") %>%
  dplyr::mutate(CVEGEO = as.factor(CVEGEO)) %>%
  dplyr::group_by(CVEGEO) %>% 
  dplyr::summarise(count = n()) %>%
  dplyr::rename(cases = count) # number of cases by municipality (CVEGEO)

# df_mxmunicipio$value <-  df_mxmunicipio$indigenous / df_mxmunicipio$pop * 100
df_mxmunicipio_journ <- df_mxmunicipio %>%
  mutate(region = as.factor(region)) %>%
  # unir tablas
  left_join(murders_by_mun,
            # indicar expl?citamente las columnas ?ndice
            by = c("region" = "CVEGEO")) %>%
  mutate(value = ifelse(is.na(cases), 0, cases))


## Complete municipalities not included
# Define the specific region codes
specific_codes <- c("02006", "04012", "07120", "07121", "07122", "07123", 
                    "07124", "07125", "17034", "17035", "17036", "23011")

# Create a list of new rows with the specified region codes and value 0
new_rows <- lapply(specific_codes, function(code) {
  data.frame(region = code, value = 0)
})

# Add the new rows to the df_mxmunicipio_journ dataset
df_mxmunicipio_journ <- bind_rows(df_mxmunicipio_journ, new_rows)

#----------------------------------------------
# MAIN PAPER'S FIGURES
#----------------------------------------------

#### Map 1
labels <- df_mxmunicipio_journ %>%
  group_by(state_code) %>%
  # Add states we mention in the text
  filter(value > 3 | state_name =="Sinaloa" | state_name =="Quintana Roo" |
           state_name =="Michoacán" |
           state_name =="Morelos" | state_name =="Ciudad de México" |
           state_name =="México" | state_name =="Oaxaca") %>% 
  mutate(group = NA) %>%
  distinct(state_code, .keep_all = T) %>%
  mutate(state_name =
           ifelse(state_name == "Ciudad de México", "Mexico City",
                  ifelse(state_name == "México", "State of Mexico", state_name)))

personalized_nudge_y <- c(
  "Chihuahua" = -.9,        # No adjustment for Chihuahua
  "Mexico City" = -0.15,      # No adjustment for Mexico City
  "Guerrero" = 0.5,       # Adjust Guerrero to the right
  "State of Mexico" = .3,  # No adjustment for State of Mexico
  "Michoacán" = -0.8,        # No adjustment for Michoac?n
  "Morelos" = -0.2,          # No adjustment for Morelos
  "Oaxaca" = -0.5,           # No adjustment for Oaxaca
  "Quintana Roo" = -1.5,     # No adjustment for Quintana Roo
  "Sinaloa" = -1,       # Adjust Sinaloa to the left
  "Tabasco" = .1,          # No adjustment for Tabasco
  "Tamaulipas" = -2,     # Adjust Tamaulipas to the right
  "Veracruz" = 0        # Adjust Veracruz to the right
)

# Define personalized adjustments for each state
personalized_nudge_x <- c(
  "Chihuahua" = 1,        # No adjustment for Chihuahua
  "Mexico City" = 0,      # No adjustment for Mexico City
  "Guerrero" = -1.05,       # Adjust Guerrero to the right
  "State of Mexico" = 0.2,  # No adjustment for State of Mexico
  "Michoacán" = -.82,        # No adjustment for Michoac?n
  "Morelos" = 0.4,          # No adjustment for Morelos
  "Oaxaca" = 0,           # No adjustment for Oaxaca
  "Quintana Roo" = 0,     # No adjustment for Quintana Roo
  "Sinaloa" = 0,       # Adjust Sinaloa to the left
  "Tabasco" = -.2,          # No adjustment for Tabasco
  "Tamaulipas" = 0,     # Adjust Tamaulipas to the right
  "Veracruz" = 0.2        # Adjust Veracruz to the right
)

mxmunicipio_choropleth(df_mxmunicipio_journ, 
                       num_colors = 1,
                       #title = "Journalists murdered in Mexico (1994-2019)",
                       legend = "Cases") +
  theme(plot.title = element_text(hjust = 0.5, size=40, face="bold"),
        title = element_text(size=35, face="bold"),
        #legend.title = element_text(color = "blue", size = 14),
        legend.text = element_text(size = 15)) +
  # For states' names
  geom_text_repel(data = labels, 
                  aes(long, lat, label = state_name),
                  size = 6, min.segment.length = 55,
                  fontface = "bold",
                  # direction = "x",
                  #  segment.size = 0.2,
                  box.padding = 0.5, max.overlaps = Inf,
                  force = 0,
                  nudge_x = personalized_nudge_x[labels$state_name],  # Personalize horizontal position
                  nudge_y = personalized_nudge_y[labels$state_name])  # Adjust vertical position

# save as PNG
ggsave("Map1.png", dpi = 600,
       height = 12, width = 10 * 2, bg = "white")

#### Figure 3
## Calculate if journalists are killed in largest municipalities or not
df_mxmunicipio_journ <- df_mxmunicipio_journ %>%
  group_by(municipio_code) %>%
  arrange(desc(pop)) %>%
  mutate(mun_order_pop = row_number()) %>%
  ungroup()



# Ensure mun_order_pop is a factor and collapse levels
df_mxmunicipio_journ$mun_order_pop <- as.factor(df_mxmunicipio_journ$mun_order_pop)

df_mxmunicipio_journ$mun_order_pop_fct <- fct_collapse(
  df_mxmunicipio_journ$mun_order_pop,
  ">=10" = levels(df_mxmunicipio_journ$mun_order_pop)[as.numeric(levels(df_mxmunicipio_journ$mun_order_pop)) >= 10]
)

# Plotting
df_mxmunicipio_journ %>%
  count(mun_order_pop_fct) %>%
  mutate(pct = prop.table(n)) %>%
  ggplot(aes(x = mun_order_pop_fct, y = pct, label = scales::percent(pct))) + 
  geom_col(position = 'dodge') + 
  geom_text(position = position_dodge(width = .9),    # move to center of bars
            vjust = -0.5,    # nudge above top of bar
            size = 10) +
  scale_y_continuous(labels = scales::percent,
                     limits = c(0,.35)) +
  scale_x_discrete(labels = c("1" = "1\nMost populated\nmunicipality\n", 
                              ">=10" = "≥10\nLeast populated\nmunicipalities")) +
  labs(title = "",
       x = "Municipalities Ranked by Population Within Each State",
       y = "Number of Journalists Killed") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 25, face = "bold"),
        axis.title.x = element_text(size = 44, face = "bold"),
        axis.title.y = element_text(size = 44, face = "bold"),
        axis.title = element_text(size = 35, face = "bold"),
        axis.text.x = element_text(size = 30),
        axis.text.y = element_text(size = 40),
        axis.text.y.left = element_text(size = 30, face = "bold",
                                        margin = margin(l = 20)),
        title = element_text(size = 15, face = "bold"),
        legend.position = "none")
ggsave("Figure3.jpg", 
       height = 12, width = 10 * 2.5)


#----------------------------------------------
# APPENDIX
#----------------------------------------------

#################### MAPS BY TIME PERIOD ########################
### 2001-2006
murders_by_mun01.06 <- journ %>%
  # filter time period
  dplyr::filter(year >= 2001 & year <=2006) %>%
  # include leading zeros in variable munic_id
  dplyr::mutate(munic_id = str_pad(mun_id, 3, pad = "0"),
                st_id = str_pad(state_id, 2, pad = "0")) %>%
  # concatenate state & munic IDs in a in a single variable
  unite(CVEGEO, c("st_id", "munic_id"), sep = "") %>%
  dplyr::mutate(CVEGEO = as.factor(CVEGEO)) %>%
  dplyr::group_by(CVEGEO) %>% 
  dplyr::summarise(count = n()) %>%
  dplyr::rename(cases = count) # number of cases by municipality (CVEGEO)

# df_mxmunicipio$value <-  df_mxmunicipio$indigenous / df_mxmunicipio$pop * 100
df_mxmunicipio_journ01.06 <- df_mxmunicipio %>%
  mutate(region = as.factor(region)) %>%
  # unir tablas
  left_join(murders_by_mun01.06,
            # indicar expl?citamente las columnas ?ndice
            by = c("region" = "CVEGEO")) %>%
  mutate(value = ifelse(is.na(cases), 0, cases))

# Add the new rows to the df_mxmunicipio_journ dataset
df_mxmunicipio_journ01.06 <- bind_rows(df_mxmunicipio_journ01.06, new_rows)

# Extract max number of journalists murdered during the time period
max.01.06 <- max(df_mxmunicipio_journ01.06$value)

labels01.06 <- df_mxmunicipio_journ01.06 %>%
  group_by(state_code) %>%
  filter(value >=  max.01.06 - 1) %>%
  mutate(group = NA) %>%
  distinct(state_code, .keep_all = T)

png("Figure A.V.c1.png",
    width = 465, height = 225, units='mm', res = 600)
map.01.06 <- mxmunicipio_choropleth(df_mxmunicipio_journ01.06, 
                                    num_colors = 1,
                                    legend = "Cases") +
  theme(plot.title = element_text(hjust = 0.5, size=40,face="bold"),
        title =element_text(size=35,face="bold"),
        legend.text = element_text(size = 15)) +
  # For states' names
  geom_text_repel(data = labels01.06, 
                  aes(long, lat, label = state_name),
                  size = 6,
                  fontface = "bold",
                  direction = "x",
                  segment.size = 0.2)
dev.off()

### 2007-2012
murders_by_mun07.12 <- journ %>%
  # filter time period
  dplyr::filter(year >= 2007 & year <=2012) %>%
  # include leading zeros in variable munic_id
  dplyr::mutate(munic_id = str_pad(mun_id, 3, pad = "0"),
                st_id = str_pad(state_id, 2, pad = "0")) %>%
  # concatenate state & munic IDs in a in a single variable
  unite(CVEGEO, c("st_id", "munic_id"), sep = "") %>%
  dplyr::mutate(CVEGEO = as.factor(CVEGEO)) %>%
  dplyr::group_by(CVEGEO) %>% 
  dplyr::summarise(count = n()) %>%
  dplyr::rename(cases = count) # number of cases by municipality (CVEGEO)

df_mxmunicipio_journ07.12 <- df_mxmunicipio %>%
  mutate(region = as.factor(region)) %>%
  # unir tablas
  left_join(murders_by_mun07.12,
            by = c("region" = "CVEGEO")) %>%
  mutate(value = ifelse(is.na(cases), 0, cases))

# Add the new rows to the df_mxmunicipio_journ dataset
df_mxmunicipio_journ07.12 <- bind_rows(df_mxmunicipio_journ07.12, new_rows)

# Extract max number of journalists murdered during the time period
max.07.12 <- max(df_mxmunicipio_journ07.12$value)

labels07.12<- df_mxmunicipio_journ07.12 %>%
  group_by(state_code) %>%
  filter(value >=  max.07.12 - 1) %>%
  mutate(group = NA) %>%
  distinct(state_code, .keep_all = T)

# Save map
png("Figure A.V.c2.png",
    width = 465, height = 225, units='mm', res = 600)
map.07.12 <- mxmunicipio_choropleth(df_mxmunicipio_journ07.12, 
                                    num_colors = 1,
                                    legend = "Cases") +
  theme(plot.title = element_text(hjust = 0.5, size=40,face="bold"),
        title =element_text(size=35,face="bold"),
        #legend.title = element_text(color = "blue", size = 14),
        legend.text = element_text(size = 15)) +
  # For states' names
  geom_text_repel(data = labels07.12, 
                  aes(long, lat, label = state_name),
                  size = 6, min.segment.length = 55,
                  box.padding = 0.5, max.overlaps = Inf,
                  force = 0,
                  fontface = "bold",
                  nudge_x = c("Veracruz" = .5),
                  nudge_y = c("Chihuahua" = -1.5,
                              "Veracruz" = .5)); map.07.12
dev.off()  



### 2013-2018
murders_by_mun13.18 <- journ %>%
  # filter time period
  dplyr::filter(year >= 2013 & year <=2018) %>%
  # include leading zeros in variable munic_id
  dplyr::mutate(munic_id = str_pad(mun_id, 3, pad = "0"),
                st_id = str_pad(state_id, 2, pad = "0")) %>%
  # concatenate state & munic IDs in a in a single variable
  unite(CVEGEO, c("st_id", "munic_id"), sep = "") %>%
  dplyr::mutate(CVEGEO = as.factor(CVEGEO)) %>%
  dplyr::group_by(CVEGEO) %>% 
  dplyr::summarise(count = n()) %>%
  dplyr::rename(cases = count) # number of cases by municipality (CVEGEO)

# df_mxmunicipio$value <-  df_mxmunicipio$indigenous / df_mxmunicipio$pop * 100
df_mxmunicipio_journ13.18 <- df_mxmunicipio %>%
  mutate(region = as.factor(region)) %>%
  # unir tablas
  left_join(murders_by_mun13.18,
            # indicar expl?citamente las columnas ?ndice
            by = c("region" = "CVEGEO")) %>%
  mutate(value = ifelse(is.na(cases), 0, cases))

# Add the new rows to the df_mxmunicipio_journ dataset
df_mxmunicipio_journ13.18 <- bind_rows(df_mxmunicipio_journ13.18, new_rows)

# Extract max number of journalists murdered during the time period
max.13.18 <- max(df_mxmunicipio_journ13.18$value)

labels13.18 <- df_mxmunicipio_journ13.18 %>%
  group_by(state_code) %>%
  filter(value >=  max.13.18 - 1) %>%
  mutate(group = NA) %>%
  distinct(state_code, .keep_all = T)


personalized_nudge_y <- c(
  "Chihuahua" = 2.1,        # No adjustment for Chihuahua
  "Guerrero" = 0.5,       # Adjust Guerrero to the right
  "Oaxaca" = -2,           # No adjustment for Oaxaca
  "Tabasco" = .2,          # No adjustment for Tabasco
  "Tamaulipas" = 0,     # Adjust Tamaulipas to the right
  "Veracruz" = 2        # Adjust Veracruz to the right
)

# Define personalized adjustments for each state
personalized_nudge_x <- c(
  "Chihuahua" = 1,        # No adjustment for Chihuahua
  "Guerrero" = -1.05,       # Adjust Guerrero to the right
  "Oaxaca" = 1,           # No adjustment for Oaxaca
  "Tabasco" = -.2,          # No adjustment for Tabasco
  "Tamaulipas" = 2.5,     # Adjust Tamaulipas to the right
  "Veracruz" = -.5        # Adjust Veracruz to the right
)


# Save map
png("Figure A.V.c3.png",
    width = 465, height = 225, units='mm', res = 600)
mxmunicipio_choropleth(df_mxmunicipio_journ13.18, 
                       num_colors = 1,
                       #title = "Journalists murdered in Mexico (1994-2019)",
                       legend = "Cases") +
  theme(plot.title = element_text(hjust = 0.5, size=40, face="bold"),
        title = element_text(size=35, face="bold"),
        #legend.title = element_text(color = "blue", size = 14),
        legend.text = element_text(size = 15)) +
  # For states' names
  geom_text_repel(data = labels13.18, 
                  aes(long, lat, label = state_name),
                  size = 6, min.segment.length = 55,
                  fontface = "bold",
                  force_pull = 0,
                  direction = "x",
                  #  segment.size = 0.2,
                  box.padding = 0, max.overlaps = Inf,
                  force = 0, nudge_y = personalized_nudge_y,
                  nudge_x = personalized_nudge_x)
dev.off()
